#ifndef __IMAGECORR_H__
#define __IMAGECORR_H__

#include <vector>
#include <list>
#include <utility>
#include <complex>
#include <map>

#include "itkFFTWRealToComplexConjugateImageFilter.h"
#include "itkFFTWComplexConjugateToRealImageFilter.h"
#include "itkVnlFFTRealToComplexConjugateImageFilter.h"
#include "itkVnlFFTComplexConjugateToRealImageFilter.h"
#include "vnl/vnl_matrix.h"

#include "alltypedefs.h"


using namespace std;

typedef itk::VnlFFTRealToComplexConjugateImageFilter 
                                        <float, 2> FFTFilterType;
typedef itk::VnlFFTComplexConjugateToRealImageFilter 
                                        < float, 2 > invFFTFilterType;
typedef FFTFilterType::OutputImageType ComplexImageType;
typedef itk::ImageDuplicator<ComplexImageType> ComplexImageDuplicator;

typedef itk::ImageRegionIterator< ComplexImageType > ComplexImageIterator;
typedef itk::ImageRegionIterator< ScalarImageType > ScalarImageIterator;

//
//
//  Region
//
//

class ImageCorrelation 
{
   public:
      enum FillType { BLACK, WHITE, NOISE, FLOOD };

   private:

      ScalarImageType::Pointer p_image1;
      ScalarImageType::Pointer p_image2;
      ScalarImageType::Pointer p_unnorm_corr_image;
      ScalarImageType::Pointer p_norm_corr_image;

      void pad( FillType fillType ) ;

      ComplexImageType::Pointer calcComplexConjugate
                   (const ComplexImageType::Pointer& I) const;

      void calcNormalizedCorrelation
         (const ComplexImageType::Pointer& F1,
          const ComplexImageType::Pointer& ccF1,
          const ComplexImageType::Pointer& F2,
          const ComplexImageType::Pointer& ccF2,
          ComplexImageType::Pointer& corr_image,
          ComplexImageType::Pointer& norm_corr_image,
          const float& scalar_norm ) const;

      ComplexImageType::Pointer calcSqrt
         (const ComplexImageType::Pointer& I) const;
      ComplexImageType::Pointer calcDiv
         (const ComplexImageType::Pointer& top, 
          const ComplexImageType::Pointer& bottom) const;
      ComplexImageType::Pointer calcMult
         (const ComplexImageType::Pointer& top, 
          const ComplexImageType::Pointer& bottom) const;
      ComplexImageType::Pointer calcDiv
         (const ComplexImageType::Pointer& top, 
          const ScalarImageType::Pointer& bottom) const;
      ComplexImageType::Pointer calcDiv
         ( const ComplexImageType::Pointer& top, 
           const float& val ) const;

      ScalarImageType::Pointer calcDenominator
         (const ComplexImageType::Pointer& F1,
          const ComplexImageType::Pointer& ccF1,
          const ComplexImageType::Pointer& F2,
          const ComplexImageType::Pointer& ccF2) const;

      void printComplexImageReal
                  (const ComplexImageType::Pointer I,
                   const long int& label) const;
      void printComplexImageImaginary
                  (const ComplexImageType::Pointer I,
                   const long int& label) const;
      void printFFT( const ComplexImageType::Pointer& I,
            const long int& label, string s = "" ) const;
      void printFFTPhase( ComplexImageType::Pointer I,
            const long int& label, string s = "" ) const;
      float calcMax(const ScalarImageType::Pointer img);
      float calcMean(const ScalarImageType::Pointer img);
      void normalize(ComplexImageType::Pointer& img) const;
      void normalize(ScalarImageType::Pointer& img) const;
      float norm( const ScalarImageType::Pointer& img ) const;

   public:


      ImageCorrelation( const ScalarImageType::Pointer& I1,
                        const ScalarImageType::Pointer& I2,
                        FillType fillType = BLACK ); 

      float getUnnormalizedMaxCorrelation();
      float getNormalizedMaxCorrelation();
};


#endif
